DNA: Prom rigid base— pairs to semiflexible polymers 
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The sequence-dependent elasticity of double-helical DNA on a nm length scale can be captured by 
the rigid base-pair model, whose strains are the relative position and orientation of adjacent base- 
pairs. Corresponding elastic potentials have been obtained from all-atom MD simulation and from 
high-resolution structural data. On the scale of a hundred nm, DNA is successfully described by a 
continuous worm-like chain model with homogeneous elastic properties characterized by a set of four 
elastic constants, which have been directly measured in single-molecule experiments. We present 
here a theory that links these experiments on different scales, by systematically coarse-graining the 
rigid base-pair model for random sequence DNA to an effective worm-like chain description. The 
average helical geometry of the molecule is exactly taken into account in our approach. We find 
that the available microscopic parameters sets predict qualitatively similar mesoscopic parameters. 
The thermal bending and twisting persistence lengths computed from MD data are 42 and 48 
nm, respectively. The static persistence lengths are generally much higher, in agreement with 
cyclization experiments. All microscopic parameter sets predict negative twist-stretch coupling. 
The variability and anisotropy of bending stiffness in short random chains lead to non-Gaussian 
bend angle distributions, but become unimportant after two helical turns. 



I. INTRODUCTION 



The sequence-dependent clastic properties of the DNA 
play a vital role in basic biological processes such as chro- 
matin organization fil 0] and gene regulation, via indi- 
rect readout [1,0,1,101 or via DNA looping 0,11 @|. The 
structure and elasticity of double helical DNA on the nm- 
scale is often described using rigid base-pair chain (RBC) 
models, in which the relative orientation and translation 
of adjacent base-pairs (bp) specify the conformation of 
the molecule EH ■ Parameter sets for rigid base-pair 
step elastic potentials were obtained from molecular dy- 
namics simulation [l2| and from an analysis of high res- 
olution crystal structure data [l3[ . We have found quali- 
tative but not quantitative agreement between these dif- 
ferent potentials in a recent study on indirect readout in 
protein-DNA binding [l4|- 

On a mesoscopic length scale, it is possible to directly 
measure force-extension relations for DNA in single- 
molecule experiments [l5|. For small external forces, 
DNA behaves as a worm-like chain (WLC) [l6j], i.e. an 
inextensible semiflexible polymer with a single parame- 
ter, the bending persistence length, and no explicit se- 
quence dependence. An extension of the classical WLC 
model, reflecting the chiral symmetry of the DNA double 
helix, includes coupled twisting and stretching degrees of 
freedom [I?], [H, [H H(| ■ These become important in a 
force regime where the DNA molecule is already pulled 
straight but not yet overstretched [2l[. Interestingly, re- 
cent measurements indicate that DNA overtwisting when 



stretched in the linear response regime [22|, |23( • 

In this article we establish a relation between these 
different levels of detail. Specifically, we coarse-grain a 
RBC to the WLC scale, while taking the average helical 
geometry of the chain exactly into account. As a result, 
we obtain the average helical parameters and the full set 
of stiffnesses for bend, twist, stretch, as well as twist- 
stretch coupling. 

It has been pointed out [24[ that the total apparent 
persistence length of a WLC is composed of a static 
part which originates from the sequence-dependent equi- 
librium bends of the molecule, and a dynamic part in- 
duced by thermal fluctuations. Their relative contribu- 
tions have been measured [H, [26[. In analogy to this 
approach, we consider the variability of static conforma- 
tions of a random RBC and from these derive the static 
and thermal persistence lengths. 

We compare the mesoscopic predictions for DNA stiff- 
ness resulting from different microscopic parametriza- 
tions in some detail, relating them to recent measure- 
ments in single-molecule experiments. 



II. RIGID BASE PAIR MODEL OF DNA 

In canonical double-stranded DNA, Watson-Crick 
base pairs are stacked into a helical column. We can 
fix a Cartesian coordinate frame to the center of each 
base pair in a standard way [13, HH , effectively averaging 
out internal distortions within the base pair. By conven- 
tion, the z-axis of this right handed orthonormal frame is 
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normal to the base pair plane and points towards the 3' 
direction of the preferred strand, while the x-axis points 
towards the major groove. 

The configuration in space of the chain is specified by 
the sequence of these frames, i.e. by a 3 x 3 rotation ma- 
trix R together with three Cartesian coordinates of the 
origin p, for each base pair step. Only for homogeneous, 
idealized and non-fluctuating B-DNA do all frames lie 
on a straight line, with their body z-axes pointing into 
a single direction. Gcncrically, the frames are displaced 
and rotated away from this idealized arrangement, due to 
both thermal fluctuations and sequence-dependent vari- 
ations in the equilibrium conformations. 

We represent the rotation and translation of the k + 1- 
th base pair frame relative to the k-th frame by a 4 x 4 
matrix, written in block form as 



9 k k+l 



Rk k+l 





Pkk+1 
1 



(1) 



Throughout the article, matrices in square brackets will 
have exactly this block structure. In idealized B-DNA 
along the z-axis, pkk+i °c cfa = (0,0,1), and Rkk+i is a 
rotation about ^3. 

This so-called homogeneous representation (see e.g. 
[29[ | has the advantage that the translation and rotation 
relating frames k and I > k can be obtained by matrix 
multiplication along the chain, 



9k I — 9kk+l9k+lk+2 ■ ' ' 9l-ll- 



(2) 



For convenience we fix the lab frame on the first base pair, 
so gik represents the frame k relative to the lab. Observe 
that Sfcfc+i = 9xk9\k+\ and g kk = e, the identity matrix. 



III. FLUCTUATIONS 

At finite temperature, a base pair step g = gkk+i in 
a RBC fluctuates around a mean or equilibrium value 
go. To parametrize these fluctuations, we first introduce 
coordinates suitable to describe small deviations from 
go- Wc will then characterize thermal fluctuations and 
the sequence randomness in terms of their second mo- 
ments. In our model, we neglect possible couplings be- 
tween neighboring base-pair steps [3(| Hl|. As will be 
explained below, the requirement of a meaningful base 
sequence nonetheless introduces some nearest-neighbor 
correlations in expectation values for random DNA. 



A. Exponential coordinates 

Any continuous group can be locally parametrized by 
its infinitesimal generators via the exponential map. In 
the c/-matrix representation, this is the ordinary matrix 
exponential exp, and the group generators {Xi} are 4x4 




FIG. 1: Frame geometry. A base pair step, connecting the 
base-pair fixed material frames e and g (left hand side). The 
frame origin trace of the corresponding screw motion is shown 
in blue. It has initial tangent v. By right multiplication with 
g ax , the same step can be described using the frames e! and g' 
(red, right hand side). They lie on the helical axis and point 
into its direction u>. 



matrices. Explicitly, in block form, 

X%= ' witn ( e i)j'fc = e i« and ( 3a ) 
x i+3 = n I* , with (di)j = Sij. (3b) 



Here, e,j& and Sij are the antisymmetric and symmetric 
tensors, respectively, and 1 < i, j, k < 3. A rotation 
around the axis is generated by Xi while a translation 
along di is generated by Xi+%. The generators satisfy 
the usual commutation relations of angular and linear 
momentum. Any group element g can be written as 
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R(0 p(0 
1 



expIfX, 



(4) 



which defines the £* as exponential coordinates of g [47| • 
The coordinate vector can be split up into two three- 
dimensional parts, £ = (uj, v). Both have a geometri- 
cal meaning: u> points along the rotation axis of R with 
||<x>|| equal to the total rotation angle, and v is the initial 
tangent j;\ p(s^), see fig. [TJ All of SE(3) except for a 
measure zero set is covered one-to-one by the coordinate 
range {£ E K 6 | |M| < tt}. 

We let the equilibrium conformation of a step go = 
exp[£gX;]. Instead of considering small additive fluc- 
tuations in £0, we also use exponential coordinates for 
go l g- Then, a fluctuating step is written as g = 
g exv[?Xi] - cxp^A^cxpieX,]. Writing £0 = 
(Tio, Roo, Two, Sho, Skj, Rio), one can check that the Tilt, 
Roll, Twist, Shift, Slide and Rise equilibrium values so 
defined have the correct symmetries under change of pre- 
ferred strand as required by the Cambridge convention 
[27l ]. This left-invariant formulation has the advantage 
that deformations with respect to different equilibrium 
positions are directly comparable and no distortions due 
to curvilinear coordinates occur. It is essential for our 
formalism which relates fluctuations given with respect 
to different frames (see below). Note however that this 
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definition of base-pair step parameters differs from those 
used in available software such as [H, [H| . We explain 
in appendix [A] how to convert between our exponential 
coordinates and the coordinate set used in (33l . |34| . A re- 
lated approach makes use of exponential coordinates for 
the rotation part of the frame transformation only [35j ] . 



B. On-axis transformation 

The screw motion s cxp[s£ l JT,;] joins the identity 
frame e with g as s increases from to 1, see fig. [TJ Its 
screw axis is determined by a vector from the origin of 
e to a point on the axis, given by p ax = \\u)\\ u> x v, 
and by its direction, u. It is the 'local helical axis' [HJ 
associated with the base pair step g. When concatenating 
many identical steps g one generates a RBC with frame 
origins lying on a regular helix with this axis. 

In addition to p ax we can define a matrix R ax which 
rotates e such that uj becomes its third direction vector. 
One choice is to take p ax as the second new direction. In 
combination, we then get 



(J a 



^ax Pax 




1 

















1 



(5) 

which takes e to a frame e' = eg ax = g ax sitting on the 
helix axis with its third direction pointing along it. One 
can check that <?' = gg ax also has these properties. The 
primed, on-axis frames arc 'local helical axis systems' in 
the terminology of [32 |. In the following, we reserve the 
name g ax for that frame transformation ([5]) which takes 
the equilibrium step go onto its own axis. 



C. Thermal fluctuations and sequence randomness 

Any base pair step fluctuates in a thermal environ- 
ment. In general the thermal mean value as well as 
the covariance matrix are sequence-dependent. In order 
to study the large scale behavior of a random sequence 
chain, we include this variability as another, indepen- 
dent source of randomness in addition to the thermal 
fluctuations [24|. I.e. we consider a random sequence 
step g = go exp[£ l X.j] which fluctuates around a global, 
sequence-independent equilibrium conformation go , with 
a covariance matrix C y = resulting from both se- 

quence and thermal fluctuations. The corresponding de- 
formation probability distribution is 



P (OdVz= P (OA(Ocie---de 



(6) 



Here, p is the probability density function (pdf) and 
dV^ = A(£)d 6 (, is the invariant volume element on the 
group, which is the Jacobian factor corresponding to our 
choice of curvilinear coordinates [36[. We can approxi- 
mate A as a constant, see appendix [Bj We now calculate 
go and C from the thermal and sequence statistics. 



We first determine go such that the expectation over 
thermal and sequence randomness, (£) = 0. This is al- 
ways possible for not too wide step distributions [37j |. 
and can be implemented by a gradient search with no 
numerical problems. 

Within a regime of linear response, the deformation 
energy of a step with fixed sequence a is a quadratic func- 
tion of the deviation from the thermal equilibrium value 
(£|cr) [48[, irrespective of the detailed nature of backbone 
connections and stacking interactions. The associated 
thermal covariance matrix is 



(7) 



On the other hand, the covariance of the sequence- 
dependent thermal mean values is given by 



(8) 



where the outermost expectation is effectively taken with 
respect to the step sequence distribution p(a) . 

Since the two sources of randomness arc independent, 
their covariances add up. One computes 



c ij =( (£ - my it - my > + < mnsw)' > 

={cv) + ct 



(9) 



Given the covariance (or stiffness) matrices and equilib- 
rium values of all sixteen dinucleotide steps, and a dis- 
tribution of relative step frequencies p(a), by computing 
go and C we have characterized a thermally fluctuating 
random sequence step in terms of its center and second 
moment. 



IV. COARSE— GRAINING 

Up to this point, step deformations and therefore also 
the covariance matrices were given with respect to a ref- 
erence frame fixed to the equilibrium base-pair frame go , 
which in general is offset and tilted relative to its own 
local helical axis. To relate the RBC deformations to a 
coarse-grained WLC model, we are much more interested 
in the elastic properties of the centerline of the chain. 
Such a centerline can be taken as the local helical axis 
for every base pair step, cf. fig. Q] This has the disadvan- 
tage that for fluctuating steps, the centerline pieces do 
not form a continuous curve. On the other hand, one can 
fit a continuous centerline globally to a stretch of a RBC 
[32| . In such an approach, the centerline depends non- 
locally on the base pair step conformations, introducing 
artificial correlations on the length scale over which the 
fitting procedure extends. 

We circumvent these problems in three steps. First we 
transform all rigid base-pairs of the chain to new frames 
of reference. These are chosen such that without fluctua- 
tions, all new bp frames lie exactly on, and point in the 
direction of a single straight helical axis. We can then 
identify and average over the unwanted shear degrees of 
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freedom. In a last step, this reduced model is averaged 
over the helical phase angle and mapped to the WLC 
models. 



A. On-axis RBC 

We would like to transform small deviations from an 
equilibrium conformation go into small deviations from 
a version of go which is on-axis. Consider first a regular 
helix composed of identical go steps. As explained in 
section lTlI B[ the on-axis step between the fc-th and k+1- 
th on-axis frames is 

.90,, = {9o9axy 1 9o + 1 9ax = 9ax90gax, (10) 

where g ax is determined entirely by go, see ([5]). Since go u 
is a transformation between on-axis frames, its rotation 
and displacement vectors point along the axis, luqu = 
||a;ou||<23 and p „ = ||poii||<fe- 
For a fluctuating RBC we calculate, 

{gikgaxy 1 gik+i9ax = g~^gkk+\9ax = go»gax e *p[€ x i}gax, 

(11) 

where gkk+i = <?o exp^'A'j] is the off-axis fluctuating 
step. The three rightmost factors in (fTTj) clearly rep- 
resent the deviation from the on-axis equilibrium step 
<70n. We introduce some standard notation. The 6x6 
adjoint matrix Adg is defined for any g G SE(3) by 
gXig~ x = [kAg)\Xj. Explicitly, if g = (R,p), one finds 

a ^={p%rr)^ ^ 

written in 3 x 3 blocks. Pulling a similarity transforma- 
tion inside the exponential series we can then rewrite (jlip 
as 

gongax 1 exp[f Xi]9ax = 00m exp^X*] (13) 

Here the deviation from the on-axis equilibrium step £ M = 
Adg~z£,, has zero mean and covariance matrix 

= {CO = {kdg-^ k C kl {Kdg-^\. (14) 

The RBC composed of steps (| X3[) is an equivalent de- 
scription of the original chain, which we may call its on- 
axis version. Intuitively, to each fluctuating frame g\ k 
of the original chain, we rigidly connected a frame g' lk 
in such a way that the primed, on-axis chain fluctuates 
about a straight, but still twisted, equilibrium confor- 
mation. This is illustrated in fig. [5] The sequence- 
dependent equilibrium conformations produce an irreg- 
ular helix. Thermal fluctuations increase irregularity. 
However, when averaging over thermal and sequence fluc- 
tuations, the on-axis configuration is exactly lined up on 
a straight helical axis. Note that we had no need to com- 
pute a fluctuating axis explicitly 



seq thermal— seq 




FIG. 2: Equivalent descriptions of a particular random se- 
quence RBC. "seq": Colored blocks represent base pairs in 
their equilibrium conformations. Wireframe blocks represent 
their on-axis counterparts, "thermal+seq" : The same, but 
with added thermal fluctuations. The top views show the re- 
duced helix axis offsets of the on-axis frames. (MD parameter 
set, base pair size scaled down by 40 % for clarity, sequence 
GCGTTGFGGGCT.) 



B. Averaging over shear variables 

The on-axis RBC has the nice property that the trans- 
lational fluctuations (£|f,£,^) = i v »: v ?\) are now exactly 
transversal to the equilibrium helix axis. They are pure 
shear modes and do not contribute to compression fluc- 
tuations along the chain. Let r\ = (u^vf,) be the vector 
of the four remaining variables. Noting that the volume 
element A depends only the angular part (see appendix 
|B|). we write 

(rfrf) = [ dVAMfo,) / dttfdttfpfc.) VW, (15) 



dv„ p{ n ) 

from which one can see that the 4x4 covariance matrix 
C y = {rfrf) is the same as C„ with its v*,v^ rows and 
columns deleted. Thus, rj has a centered distribution 
with covariance matrix C . Here and in the following, 
~ indicates deletion of the shear^rows and columns in 
an on-axis, 6x6 matrix. E.g, Ad is the 4x4 adjoint 
matrix. Also unless noted otherwise, we consider only 
on-axis quantities and suppress the •„ subscript in the 
following. 

C. Correlations induced by sequence 

While we assume thermal fluctuations of neighboring 
steps to be independent random variables, there are nev- 
ertheless correlations in the sequence identity of neighbor- 
ing base-pair steps. Any realization of a random sequence 
of dinuclcotidc steps must be 'continuous', e.g. an = AG 
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implies that 023 can only start with a G. These correla- 
tions need to be taken into account when calculating ex- 
pectation values for random sequences of DNA. For this 
purpose, we now consider the combined fluctuations of a 
short RBC consisting of m base pair steps. 

Assuming independent, identically distributed bases 
we obtain the joint pdf of sequence steps along the chain 
as the product of the base pdfs. p(ai2, ■ ■ ■ , Cm-m+i) = 
rifeS. P(bk)- This implies that the covariance between 
thermal mean values, 



( (vlk+i\<Tkk+i)(vit +1 Wu+i) ) = 



(16) 

Here we introduced a nearest-neighbor term C 1 which 
will be computed below. No nearest-neighbor correlation 
occurs in the thermal covariances by assumption. 

We are now in a position to combine the m base pair 
steps of our chain into one compound step. Here, m 
must be small enough that the typical deviation angles 
of the compound step from equilibrium stay small. I.c, 
the short chain must be well approximated by a (helical) 
rigid rod. Successively commuting the on-axis equilib- 
rium steps <7on to the left and introducing the adjoint 
matrix as in (1131). one arrives at 






I = k 


< C i 


l = k + l 


Cf 


I = k- 1 


,0 


otherwise. 



Son exp[£i 2ll Xi] • • • 0on exp[£ 



m m+l 



,x t ] = 



5o,i 



5>d 5o - m+fc 6 fe+ll 



k=l 



+ O(0 2 . (17) 



The sum in parentheses is the deviation £i m +i M of the 
fluctuating compound step from its equilibrium value g^ r 

The corresponding reduced Ad matrix has a simple 
form. Using (|12[) and noting that pq u oc u>q u oc c/3 we 
obtain 



r]i m+l 



= J2 Ad 



(Jo 



r]kk+i , where 



k=l 



cos 1 1 oj 1 1 sin||wo|| 0\ 
I -sin||w || cos||w || 




(18a) 



(18b) 







1/ 



One sees that the w^' 2 components are successively ro- 
tated around the c?3 axis, while the w 3 ,v 3 components 
are unaffected. 

What is the covariance matrix Cf m+1 = 

(Vim+iVim+i) °f the compound deviation? Using 
(|16p . we are left with a sum of appropriately trans- 
formed single-step covariances C = {C a ) + C and in 
addition a sum of nearest neighbor cross-terms involving 



Ci: 



m — 1 



Cim+i = J2 Ad 5o -'CAd T 



<h 



1=0 

m-2 



J2 Ad.g 'C x AdV, 



(19) 



where C x = C 1 Ad T g' 1 + Ad g~ l C [ 



-1 *~iT 



The cross-covariance C x represents the fact that near- 
est neighbor equilibrium steps are correlated and their 
frames of reference are rotated by an angle ||<j->oII- 

Note that two neighboring compound steps are still 
correlated by sequence continuity at their interface. From 
we have the recursion relation 



Ci 



m+l 



Ad 9 o 1 C lm Ad r go 1 + C + C : 



(20) 



The same relation is obeyed by a sequence of independent 
steps with covariance matrix C = C +C x . We conclude 
that except for a boundary term C x at the beginning 
of the chain, a RBC with independent steps and with 
covariances C exhibits the same effective covariance as 
the original, short range correlated chain with C and 
C x . The relative error in effective compound covariance 
is of order 1/m. 



D. Averaging over the helical phase 

A shear-averaged, on-axis RBC still has a finite equi- 
librium twist and anisotropic bending stiffness. To relate 
it to a WLC with isotropic bending rigidity, we perform 
an average over a continuous helical phase angle rotation 
of the reference frame [38[ . An on-axis covariance matrix 
which is rotated by a helical phase angle <fi around the 
average local helical axis (see (HU), is 



C{<j>) = Ad ff0 CAd T «7 , 



(21) 



where g^ = cxplc/jX^] is a pure rotation by an angle (j> 
around d%. Since Adg^ has the form (|18b[) . the helical 
phase average is seen to be 



C 



1 

27 



c {cf>)d<i) 



C LL +C 
2 







c 



0_ 






a 34 



\ 


a 34 

c 44 

(22) 



V 



From C one can read off the bend and twist persistence 
lengths as If, = h u /C u and l t = h„/C 33 , respectively, 
where the on-axis helical rise is h n = ||pon||- The WLC 
stiffness matrix j3S = C ' ~ 1 can be found by inversion 
and has the same block structure, see also appendix [B] 
Its nonzero components are the bend, twist, stretch and 
twist-stretch coupling stiffness coefficients. 
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E. Coarse— graining relations 

We have derived all WLC elastic parameters starting 
from an arbitrarily oriented and offset RBC. We now dis- 
cuss in some detail how these coarse-grained parameters 
are related to the microscopic RBC parameters. 



1. Equilibrium step 

The transformation of the equilibrium step onto the 
helical axis (fTU|) leaves the total rotation angle invari- 
ant. Therefore the equilibrium twist of go u is 9 U = 
ll^onll = ll w oll > l w ol- I- e i the twist per base pair of 
the WLC equals the total angle of rotation, not the Tw 
angle of the off-axis step. The equilibrium rise on axis 
is h„ = UpoiiII = WqPo/II^oII which is different from both 
off-axis quantities ||po|| and Pq. These differences are of 
order 0(u>q + w 2 ) 2 so they become important only when 
the equilibrium rotation axis too has significant roll and 
tilt with respect to the material frame, i.e. when the lo- 
cal helical parameters Inclination and Tip [27| arc not 
negligible. 



2. Fluctuations 



In this expression, the rotational block Cj, is merely 
a rotated version of the off-axis rotational block C^^\ 
In contrast, the translational block C^"^ and the cou- 
pling block C^° v ^ have 'leverage terms', since rotational 
fluctuations about directions perpendicular to the offset 
vector contribute through a cross product with p ax . For 
C,, , these involve the off-axis coupling C^ VLU ^ in first 
order and rotational fluctuations in second order 

in ||Pox||- The coupling block C,, ^ has contributions 
from in first order. These leverage terms persist 

in the reduced WLC covariance matrix C . They are the 
remainder of the microscopic description of fluctuations 
with respect to a material frame that is offset from the 
average helical axis. 

Consider for example a base pair step that exhibits x- 
displacement but no Inclination or Tip, i.e. p ax oc di,co oc 
dz,R a x = h- Then ([23]) implies that any coupled Roll- 
Rise (C 26 ) and Roll (C 22 ) fluctuations will add to the 
stretching fluctuations C,f 6 of the chain. In addition, 
the off- axis Roll-Twist fluctuation (C 23 ) contributes to 
twist-stretch coupling fluctuation on axis, C 36 . 

When Inclination or Tip are nonzero, then due to the 
additional rotation R ax also Shift and Slide fluctuations 
contribute to the resulting WLC parameters. It is there- 
fore essential to transform to an on-axis frame before 
averaging over the shear degrees of freedom. 



Unlike the equilibrium step, the covariance matrix is 
changed not only by the rotation R ax but also by the 
shift p ax onto the average local helix axis. Intuitively, 
the on-axis frame g' is rigidly connected to g, cf. figs. 
[TJ Therefore, a rotational fluctuation of g with rota- 
tion vector to' will result in an additional translational 
fluctuations of g' equal to to' x p ax . 

A familiar example of this geometrical effect is the 
stretching of an ordinary coil spring along its helix 
axis. In the wire material, this deformation corresponds 
mainly to torsion, i.e. a rotational deformation of consec- 
utive wire segments. On a larger scale, this deformation 
is levered into a translation of one coil end along the he- 
lix axis. The transformation (|14p captures exactly this 
lever arm effect, which is proportional to the total axial 
displacement ||p x|| an d so becomes relevant if the chain 
deviates from an idealized B-DNA form. 



We calculate explicitly the 3x3 blocks C,) ; of C„, (I3J), 
in terms of the corresponding blocks C i - ab ^ > of C, using 
dH]) and (fT2"|) . Here a,b £ {u>, v} stand for the set of ro- 
tational or translational components, respectively. Fur- 



ther, we let C^' = R^C^Rar, and PL 



R 



an antisymmetric matrix. Using this notation, 



ax jPax e ii 



f^{vv)t pi fi{u)U))t pi 

/~i{vuj)i pf ax L - / r ax 

i ax ^ \(~i(vuj)tpt pi /~i{ljv)i 



(23) 



F. Numerical verification 

We tested our coarse-graining relations by performing 
a simple-sampling Monte Carlo simulation. After gener- 
ating a random sequence, for each dinucleotide, random 
conformations were drawn according to a Gaussian distri- 
bution with the corresponding microscopic parameters. 
The measured mean squared base-pair center end-to- 
end distances are shown in fig. [3J The theoretical curves 
(R 2 ) = 2ll b -2ll(l-e- l / lb ) for an inextensible WLC using 
the computed contour and bending persistence lengths, / 
and lb, fit the simulation data to within numerical error. 
The only exceptions occur below 3 nm, where the inex- 
tensible WLC model is not a good description for the full 
shear able helical RBC. 



V. WLC PARAMETERS OF DIFFERENT RBC 
PARAMETER SETS 

As a result of the coarse-graining procedure outlined 
above, we obtain a set of WLC parameters from a set of 
sequence-dependent RBC stiffness (or covariance) ma- 
trices and equilibrium offsets. There are several different 
parameter sets available in the literature, extracted from 
analysis of X-ray crystal structures of DNA [l3| and from 
molecular dynamics simulation [l2l . |39| . 

For the stiffnesses obtained from structural data, the 
missing thermal energy scale is substituted by an "ef- 
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FIG. 3: Comparison of a simple sampling Monte Carlo simu- 
lation of random-sequence DNA to our theory. Symbols des- 
ignate the measured mean squared end-to-end distances for 
static disorder only (upper row) and for static plus thermal 
fluctuations (lower row). The theoretical curves assuming a 
WLC model for static disorder, uncorrelated static disorder, 
and for static plus thermal fluctuations are shown in blue, red 
and orange, respectively. MD microscopic parameter set, as 
explained below. 





2n 7, 


7 7 

lb H 


on 11 /Q o 33 /□ c* 44 c 34 
pj po po pj 


B 

P 
MD 
MP 


10.1 0.334 
10.5 0.334 
11.9 0.318 
10.5 0.334 


27.1 15.2 
43.4 35.7 
38.9 45.1 
42.8 47.8 


81.1 46.7 1300. -39.9 
130. 117. 1280. -116. 
122. 158. 586. -96.3 
128. 150. 1020. -81. 


units 


1 nm 


nm nm 


rad~ 2 rad -2 nm~ 2 (nmrad) -1 



TABLE I: WLC geometry, persistence lengths and stiffness 
parameters for the considered potentials. In our units, PS 11 
and PS 33 are the bending and twisting persistence lengths, 
given in base pairs. 





PS 11 PS 33 PS M PS 3A 


Gore et aZ.[23] 
Lionnet et al, [221 


163 ± 15 327 ± 15 781 ± 150 -64 ± 15 
294 710 -47 ± 20 


units 


rad -2 rad -2 nm -2 (nmrad) -1 



TABLE II: Experimental stiffness parameters as given in the 
literature, converted to our single-step units. The conversion 
factor for B, C, G, S from [U is P/h„. The conversion factors 
forB,C",Din Q are respectively, Bfjh 3 , l/h„ , (9,,/h 2 . Beware 
of a missing 1/2 factor in their first formula. 



fective temperature". We here use the effective temper- 
atures determined in a previous study [l4| by equating 
the total, microscopic fluctuation strengths of the crys- 
tal and MD covariance matrices. The absolute magni- 
tudes of all parameters derived from structural data (B 
for B-DNA crystal and P for Protcin»DNA cocrystals) 
arc therefore depend on our choice of effective temper- 
ature. Still, their relative magnitudes are properties of 
the microscopic structural data set independent of this 
choice. No such restrictions apply to the MD parameters 
(MD), since here the temperature is set by the simula- 
tion. We also include a hybrid parametrization (MP) 
which combines the equilibrium values from the P»DNA 
dataset with the stiffness matrices from MD. This com- 
bined potential compared favorably to the others in bind- 
ing affinity prediction [l4| . It can be seen as a version of 
the MD potential which is corrected for the well known 
undertwist occurring in MD simulations. For MD and 
MP, our coarse-graining involves no free parameter. 

In table U we show the resulting WLC stiffness param- 
eters and geometry. For the crystal parameter sets, the 
equilibrium rise and twist are close to the commonly ac- 
cepted values of 0.34 nm/step and 10.5 bp/turn. The 
MD rise and twist are both l ow, a known effect for the 
force field used in that study [40j . 

The MD bending persistence length is smaller than the 
commonly accepted values at physiological conditions, 
see e.g. [23]. It is also somewhat below the range of 
45 - 47 nm found experimentally [U S3, EH at the con- 
ditions of the simulation of ~ 100 mM Na + . (However, 
in [44| a lower experimental value is reported.) The low 
equilibrium Rise of the MD conformations accounts for 
half of this deviation. 



The twisting persistence lengths of all parameters sets 
are similar to the bending persistence lengths, which is 
in stark contrast to measurements of twisting persistence 
in single-molecule studies which give a value close to 100 
nm, see [ll| for a review. For the crystal parameter sets 
one might argue that this indicates that torsional defor- 
mations carry more elastic energy than bending defor- 
mations, thus 'violating' an assumed cquipartition of en- 
ergy. However, for the MD parameter set, this is clearly 
not the case; the simulated DNA oligomers were indeed 
more twistable than experimental values for DNA sug- 
gest. 

The twist-stretch coupling is negative in all cases. This 
is counter-intuitive since it implies that DNA overwinds 
in linear response to stretching. The same sign of the 
coupling is also found in the "naive" Twist-Rise coupling 
stiffness of the original stiffness matrices for (8, 9, 10) of 
the 10 unique basepair steps in the (B, P, MD) param- 
eters, respectively. Negative twist-stretch coupling has 
recently been observed in single-molecule experiments at 
low applied tension [2^, HH . We show the full elastic pa- 
rameters collected in these articles in table [H] for compar- 
ison. Generally, the agreement between the microscopic 
parameter sets and single-molecule data is better for the 
twist-stretch coupling than for the twisting rigidity. The 
stretching modulus j3S M differs by about a factor of 2 
between the crystal and MD parameters, with the exper- 
imental value inbetween. We remark that no rescaling 
by a different effective temperature can bring all crystal 
stiffness parameters into reasonable agreement with ex- 
periment since the various deviations occur in opposite 
directions. 

Instead of looking at random DNA, we can consider 
thermal and sequence fluctuations separately. Table IIIII 
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full 


lb/nm 
thermal static 


static' 


full 


It /nm 
thermal static static' 


B 


27.1 


29.5 


327. 


211. 


15.2 


15.4 


1260 


88.3 


P 


43.4 


45.3 


1040. 


575. 


35.7 


36.3 


2430 


172. 


MD 


38.9 


42. 


519. 


175. 


45.1 


47.7 


818. 


256. 


MP 


42.8 


44.6 


1040. 


575. 


47.8 


48.8 


2340. 


172. 



TABLE III: Thermal and static contributions to the apparent 
persistence length for different potentials. For comparison, 
the static' column shows the static persistence lengths when 
sequence continuity is disregarded. 





AA AC AG AT GG CG 


B 

h P 

nm MD 

MP 


32.2 18.8 43.2 35. 47.5 27.1 
40.2 49.5 49.2 40.1 48.7 40.2 
45.9 40.2 45.7 33. 51.9 38.4 
47. 44.1 46.3 37. 53.8 42.1 


B 

U P 

nm MD 

MP 


9.4 7.35 25.6 18.7 19.1 26. 

53.5 29. 44.2 34.9 32. 34.6 
46.1 42.3 48.6 61.8 59.7 38.2 

45.6 44.2 50. 63. 60.4 40.2 



TABLE IV: Comparison of persistence lengths of all six 
unique repetitive sequences of period two, for the MP and 
MD parametrizations. 



shows the corresponding static and thermal persistence 
lengths [24|. It follows from eqn. ([9]) that their inverses 
add up to give the inverse apparent (or random DNA) 
persistence length. In disagreement with the cryo-EM 
study 25j we find that the static persistence lengths are 
much higher than the thermal ones, leading to a correc- 
tion of only a few nm random DNA persistence lengths. 
This is in accordance with cyclization data [2^|. Also, 
the static lb for the P parameter sets correctly reproduces 
the value found numerically in that study, using the same 
parameter set. When we disregard the requirement of se- 
quence continuity by setting to zero all C \ contributions 
in (|16[) , static variability is strongly overestimated (more 
than tenfold for twist). 

The range over which the stiffness of random B-DNA 
can vary depending on sequence can be estimated from 
the persistence lengths of all six unique repetitive se- 
quences of period 2, given in table HVi see also (|24|) . 
Generally, 4 has similar dependence on the sequence in 
all considered potentials, while the predictions for l t are 
less correlated. The large deviations in the B-DNA pa- 
rameter set are likely due to insufficient statistics 1 1 31 ) - 
The TA(=AT) repeat stands out as the most bendable 
sequence which is at the same time torsionally stiff. An- 
other common trend is that poly-G DNA is compara- 
tively stiff with respect to bending. 

A more detailed view of the sequence variability of 
WLC stiffness is given in table [V] for the MP hybrid 
potential. The stretch modulus and the twist-stretch 
coupling depend on the sequence in a correlated way. 
The rightmost column shows the ratio of overtwist over 
elongation in response to an external stretching force, 
Tresp = C 34 /C M . When a repetitive sequence is cut 
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pS 


ps 33 


ps 44 


n ry 34 

ps 


Tresp 


AA 


144. 


141. 


976. 


-38.3 


0.27 


AC 


132. 


142. 


1140. 


-105. 


0.74 


AG 


139. 


159. 


1120. 


-103. 


0.64 


AT 


111. 


195. 


975. 


-80.1 


0.41 


GG 


159. 


186. 


1090. 


-89.9 


0.48 


CG 


124. 


126. 


831. 


-78.5 


0.62 


units 


rad -2 


rad~ 2 


nm~ 2 


(nmrad) 


rad/nm 



TABLE V: Comparison of stiffness parameters of all six 
unique repetitive sequences of period two, for the MP hybrid 
parametrization. 



by one bp and then stretched to the original length, the 
"missing twist" at the last bp ranges from 29 (AA) to 20 
(AC) degrees undertwist. 



VI. VARIABILITY OF STIFFNESS 
A. Bend angle distributions for short chains 

The combined covariance matrix C i m +i gives the 
second moment of the distribution pi m +i of deforma- 
tions, observed in a thermal ensemble of random sequence 
oligonucleotides of length m steps. Here it is not neces- 
sary that the single step deformation distributions have a 
Gaussian shape. Indeed such an assumption depends on 
the choice of coordinates, and is not justified by experi- 
ments. Nevertheless, let us for the moment additionally 
assume that the single step deformation distributions are 
in fact Gaussians. In that case the deformation of a spe- 
cific compound step again follows a Gaussian distribu- 
tion p(rji m +i |<7i m+i), since it is the result of a convolu- 
tion. However even in this case, after averaging over se- 
quence randomness, the corresponding deformation dis- 
tribution of a random compound step p\ m +i(jli m+i)j 
deviates from a Gaussian shape. This comes about by 
averaging together several Gaussians with different off- 
sets and widths. To illustrate this point, we show in 
fig. 2] the effective potential U e g for the total bend angle 
■d = ((?7i im+ i) 2 + (vf^n+i) 2 ) 1 ^ 2 01 random sequence com- 
pound steps of different lengths. It is extracted from his- 
tograms of a simulation as described in section HVFl For 
compound steps shorter that 5 bp, the effective potentials 
stay well below their respective harmonic approximations 
which are tailored to reproduce the second moment of the 
bend angle distribution. These second moments agree to 
within 1% with the bend angle variance of a WLC model 
with the same persistence length, again confirming our 
calculation. Thus for short random chains, large bend- 
ing angles occur much more frequently than would be 
expected from a WLC model with matching persistence 
length. This effect is the combined result of the varying 
bending stiffness coming from sequence as well as from 
anisotropic bending (see below). For the parametriza- 
tions we considered, this effect is very small for compound 



0. 0.05TT 0.1TT 0.15TT 0.2TT 0.25TT 0.3TT 

$/rad 

FIG. 4: Effective potential for the total bend angle ■& (green, 
symbols). The blue curves show the harmonic approximation 
to the effective potential that yields the same variance ($ 2 ). 
Compound step length, from left to right: 1,2,3,5,10 bp. MP 
parameter set. 



steps above 5 bp and is thus insufficient to explain the 
frequent large bending angles observed in a recent AFM 
study of DNA adsorbed on a surface [45[ , on length scales 
of more than 15 bp. 

We note that the shape of the bend angle distribution 
depends on what exactly is considered the local bend an- 
gle. Instead of $ as defined above one could take the 
angle between vectors (pi+i — Pi) connecting successive 
bp centers. For this choice, below 5 bp steps, the oppo- 
site behavior is seen: The second moment is increased 
while extreme bend angles are suppressed compared to 
the WLC prediction (data not shown), although both 
bend angle definitions agree on scales longer than a heli- 
cal turn. 



B. Decay of variability 

To investigate the shape of pi m +i for small m in some 
more detail, first consider a compound step with a fixed 
m step sequence <j\ m +i = (uu, • ■ • , ff mra +i)- Essentially 
just by disallowing sequence randomness in (|19p . we com- 
pute the combined covariance matrix of this compound 
step to be 

m — 1 

C aim+1 = ]T Ad ff0 -'C ffw . lB+ , I AdV 1 (24) 

1=0 

valid for small deviations relative to the sequence- 
averaged equilibrium step g™. 

We can describe the sequence variability of compound 
step covariances in terms of their first moments in the 
same way as was done for the equilibrium conforma- 
tions in section IIII CI While the mean covariance matrix 
A/ y = (C l J lm is just the sequence average of (IM1) . 



bp 

FIG. 5: Relative spread Al/l of the bending (Z&, green) and 
twist (h, blue) persistence lengths vs. compound step length. 
Ignoring sequence continuity leads to overestimation of the 
stiffness variability (triangles). 

the covariances of the entries of the thermal covariance 
matrix are given by 

= <(C« m+I - Af«)(C£ ra+1 - M kl )). (25) 

This expression can be split into diagonal and nearest 
neighbor terms in analogy to (|16p . again reflecting se- 
quence continuity. In particular, it is impossible to com- 
bine two of the comparatively soft pyrimidine-purine [l3j 
steps in a row. The resulting relative spread of thermal 
persistence lengths among random sequence compound 
steps is shown in fig. [5l Explicitly, Alb = (Vim+i) 1 ^ 2 
and Alt = (^i 3 m+i) 1/2 - Note that already after one full 
turn, variability in stiffness is down to 5 %, and that se- 
quence continuity results in reduced variability compared 
to a model with independent step sequences. 

In summary, we remark that the detailed shape of the 
deformation distributions is not known, and there is no 
reason to believe it should be Gaussian for small step 
numbers. Even when starting with Gaussians for the 
single steps, we obyain clearly non-Gaussian shapes for 
the random sequence bend angle distributions up to a few 
steps. For the long-wavelength behavior of the chain, the 
relevant quantities are just the first and second moments 
which we have calculated in section [TV] 

C. Anisotropic bending 

Another feature of short compound steps is their 
anisotropic bending stiffness. It is clear that on scales 
much longer than a full turn, the molecule behaves as a 
uniformly bending rod, at least for small deformations. 
Using the compound covariance C (before helical averag- 
ing is performed, see fH])) we can quantify the decay of 
anisotropy for random sequence chains. The ratio of the 
principal bending stiffnesses as a function of chain length 
is shown in fig. El Since linear response is always symmet- 
ric, bending into major and minor groove has the same 
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FIG. 6: Bending anisotropy. The ratio of larger over smaller 
bending stiffness decays in an oscillating fashion with com- 
pound step length. MP parameter set. 

stiffness for small deformations. As a result, the bending 
anisotropy has minima every half turn of the double he- 
lix. Since the 21 bp chain has exactly two full turns, the 
anisotropy is suppressed completely, but also already a 5 
bp compound step is not far from 5.25 bp and behaves 
essentially isotropic. 



VII. CONCLUSIONS 

In this article we have shown a way to quantita- 
tively connect experiments on DNA elasticity on differ- 
ent length scales. Starting from atomistic data, DNA 
deformations are described in terms of a rigid base-pair 
chain model in a first step [36| . We then relate the stiff- 
ness expressed in terms of rigid base-pair deformations, 
to the long-wavelength WLC parameters of a random 
chain. In this coarse-graining step it is essential to prop- 
erly account for the helical base-pair geometry. For this 
purpose we introduce an on-axis version of the rigid base- 
pair chain, which on average has ideal B-DNA shape. 
This makes it straightforward to integrate over the shear 
degrees of freedom and helical phase, to finally obtain all 
four linear elastic constants allowed by the large-scale 
symmetry of the molecule [H, H, ^ • 

Our results allow a direct comparison of the different 
microscopic effective potentials to single molecule and cy- 
clization experiments. It involves no free parameter for 
MD simulation data, and a single parameter (the effec- 
tive temperature) for structural data. Wc find good qual- 
itative agreement, including the negative sign of twist- 
stretch coupling. Quantitatively, the microscopic bend- 
ing persistence lengths agree best with recent single- 
molecule data. The twist persistence is about 50 % lower, 
and the magnitudes of compressional modulus and twist- 
stretch coupling are roughly 50 % higher than the meso- 
scopic experimental values. 

Does the involved computation of macroscopic param- 
eters actually make a noticeable difference? The calcula- 



te 12. -1. 4. -23. 

P 8. -7. 8. 24. 

MD 11. -9. 67. 52. 

MP 6. -4. -3. 44. 



TABLE VI: Relative error in stiffness parameters made when 
using naive matrix elements instead of the coarse-grained pa- 
rameters described above. Values are given in per cent. 

tions can be simplified in two ways: By disregarding the 
details of average helical geometry of the chain, and by 
treating the base sequence of adjacent steps as indepen- 
dent, i.e. disregarding sequence continuity. 

If the average helix geometry is treated correctly but 
sequence continuity is disregarded, static variability is 
strongly overestimated (table llllj) . Overall this remains 
a minor effect since the thermal fluctuations dominate. In 
addition, such a simplification leads to an overestimation 
of stiffness variability for short random oligonucleotides, 
see fig. [5j 

On the other hand one can exclude static variability 
and treat the helix geometry as ideal B-DNA from the 
beginning. Starting from the sequence-averaged, off-axis 
covariance matrix, one would perform an average over 
Shift, Slide and helical phase angle and invert to get an 
"naive" stiffness matrix S na . The relative error made in 
such a naive computation, e lJ = (S% a — S y ] )/S %3 is shown 
in table I VII While the bending and twisting stiffnesses 
are well approximated by the naive guess, the error in 
stretch modulus and twist-stretch coupling is consider- 
able. For these terms, leverage due to the axis offset be- 
comes important as explained in scction HVEl Especially 
the naive twist-stretch coupling is not negative enough. 

The procedure wc describe involves no approximations 
regarding the geometry. This makes it directly applicable 
to alternative DNA structures, once microscopic covari- 
ance matrices are available. In fact, the more the average 
geometry deviates from idealized B-DNA, the greater is 
the need to treat the helical geometry correctly. Already 
for the MD parameter set, the error when using the naive 
geometry is quite important. 

The main model assumption is that thermal defor- 
mation fluctuations of ncighoring steps are independent . 
Another limitation of any rigid base-pair model is that 
internal deformation fluctuations of a base-pair such as 
propeller twist or buckle, are not explicit and thus effec- 
tively treated as uncorrelated between base pairs. 

Our framework can be extended to improve on both of 
these points. Nearest-neighbor correlations in base-pair 
parameters may be included by extending the model to 
a full Markov chain. Internal deformations could then 
be added by extending the configuration space, leading 
to a bi-rod (4o| in the continuum limit. However for ei- 
ther of these interesting generalizations, a microscopic 
paramctrization is an open challenge in itself. The fact 
that dinucleotide step stiffness depends overall rather 
weakly on the flanking sequence [30| and the encouraging 
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agreement with mesoscopic data we found, suggest that 
the main features of coarse-grained DNA elasticity are 
captured already by our more basic model. 

In view of an experimental precision of the order of 
one percent for the mesocopic bending rigidity (26[, we 
consider a quantitatively correct relation between meso- 
scopic and microscopic stiffness parameters essential. We 
hope that the method presented in this article proves use- 
ful in providing this link. 
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APPENDIX A: COORDINATE CONVERSION 

How does one obtain the covariance matrix C and 
equilibrium conformations go for a given collection 
{gk}i<k<N of bp frame conformations? We can first 
determine go by requiring that {<7o~ 1 <7fc} has mean in 
exponential coordinates. For not too wide distributions, 
such a center always exists and is unique [oTI ]. Then, 
(jij — ) is the standard covariance matrix of {g ~ 1 gk} 
in exponential coordinates. 

However, for the potential parametrizations considered 
here, only the equilibrium values Co and covariance ma- 
trices C l ? = ((C~Co) l (C - CoP) with respect to the global 
coordinates C = (fl,T, p,q 1 ,q 2 ,q3) as defined in [34[ and 
used in [33j], are given. Here, 6 = (fi, t, p) are Twist, Tilt 
and Roll angles but differ from our choice of angles. The 
Q — (?i j 92, 93) gives the translation vector with respect 
to the mid- frame R m . The conversion formulas are, 



R 



R(C) = exp((fi/2 - arctan(r/p))e 3 ) exp(vV + r 2 e 2 ) 
cxp((f2/2 + arctan(r/p))e3), 

(C) =cxp((fi/2 - arctan(r/ / 9))e 3 ) cxp(vV + r 2 /2e 2 ) 
exp((arctan(r/p))e3), and 

p(C) = Rm(t)q, 

(Al) 

together determining the frame conformation 17(C)- We 
checked that the variation of the volume element in the 
region of noticeable probability around go is small com- 
pared to the variations in the probability density. There- 
fore neglecting the former, we get go = g(Co)- I n linear 
order around the equilibrium position, we can then trans- 
form the covariance matrix given in (^-coordinates to 
exponential coordinates using just the Jacobian matrix 



Jo of the coordinate transition map C 1— > C(C) = log(g(Q). 
This gives C = JoCqJq . We have calculated Jo = g||^ 
analytically. Its 3 x 3 blocks are 



=l/2tv{eiR T d ej R) 
=0 



du) 1 
W 
do/ 

-QQ - =(R d e jR m idq) 1 

d(q) 



(A2) 



All coarse graining calculations presented in this article 
use the matrices C converted in this way as a starting 
point. 

The exponential coordinates of the equilibrium confor- 
mations have the usual symmetries under strand change 
and reading direction reversal: Denote by a the se- 
quence complementary to a, e.g. AG = CT, and let 
E = diag(-l, 1,1, -1,1,1). Then as a -> a, Co = 
(Tio,Roo,Tw ,Sh ,Slo,Rio) -> Duc to thc & de- 

pendent coordinate conversion above, the body-frame 
covariance matrix does not obey the corresponding sym- 
metries, C ECE. While this may seem a serious 
drawback of the coordinate system we use here, it turns 
out that in the on-axis, shear and helical phase averaged 
covariance matrices, the strand-exchange symmetry is 
re-established. Therefore, our coarse-grained results are 
indeed independent of the reading sense. 



APPENDIX B: VOLUME ELEMENT 



In our coordinates, In /1(C) = — ■ 
that in a Gaussian approximation, 



^ll 2 + o(|M| 4 ) 



> +A >^ 3 d 6 C, A 



kh 



a 



(Bl) 

Here, I3 and O3 are the 3x3 identity and zero matrices, 
respectively. In DNA, the distributions of single 
steps are very narrow. Therefore when computing mo- 
ments, in particular the covariance matrix C 1 ^ = (C 1 C j )j 
we can extend the integration boundaries to infinity with 
negligible error. Performing the integral we then get the 
relation (3S + A = C _1 . Since j3S 3> A, in making the 
approximation {3S = C~ l , we introduce an error of less 
than 1% for typical B-DNA steps. I.e. thc stiffness matrix 
/3S is indeed given by the inverse of the covariance. 
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